      subroutine sun(y,m,DD,UT,lon,lat,RA,Dec,LST,Az,El,mjd,day)

      implicit none

      integer y                         !Year
      integer m                         !Month
      integer DD                        !Day
      integer mjd                       !Modified Julian Date
      real UT                           !UTC in hours
      real RA,Dec                       !RA and Dec of sun

C  NB: Double caps here are single caps in the writeup.

C  Orbital elements of the Sun (also N=0, i=0, a=1):
      real w                            !Argument of perihelion
      real e                            !Eccentricity
      real MM                           !Mean anomaly
      real Ls                           !Mean longitude

C  Other standard variables:
      real v                            !True anomaly
      real EE                           !Eccentric anomaly
      real ecl                          !Obliquity of the ecliptic
      real d                            !Ephemeris time argument in days
      real r                            !Distance to sun, AU
      real xv,yv                        !x and y coords in ecliptic
      real lonsun                       !Ecliptic long and lat of sun
C Ecliptic coords of sun (geocentric)
      real xs,ys
C Equatorial coords of sun (geocentric)
      real xe,ye,ze
      real lon,lat
      real GMST0,LST,HA
      real xx,yy,zz
      real xhor,yhor,zhor
      real Az,El

      real day
      real rad
      data rad/57.2957795/

C  Time in days, with Jan 0, 2000 equal to 0.0:
      d=367*y - 7*(y+(m+9)/12)/4 + 275*m/9 + DD - 730530 + UT/24.0
      mjd=d + 51543
      ecl = 23.4393 - 3.563e-7 * d

C  Compute updated orbital elements for Sun:
      w = 282.9404 + 4.70935e-5 * d
      e = 0.016709 - 1.151e-9 * d
      MM = mod(356.0470d0 + 0.9856002585d0 * d + 360000.d0,360.d0)
      Ls = mod(w+MM+720.0,360.0)

      EE = MM + e*rad*sin(MM/rad) * (1.0 + e*cos(M/rad))
      EE = EE - (EE - e*rad*sin(EE/rad)-MM) / (1.0 - e*cos(EE/rad))

      xv = cos(EE/rad) - e
      yv = sqrt(1.0-e*e) * sin(EE/rad)
      v = rad*atan2(yv,xv)
      r = sqrt(xv*xv + yv*yv)
      lonsun = mod(v + w + 720.0,360.0)
C  Ecliptic coordinates of sun (rectangular):
      xs = r * cos(lonsun/rad)
      ys = r * sin(lonsun/rad)

C  Equatorial coordinates of sun (rectangular):
      xe = xs
      ye = ys * cos(ecl/rad)
      ze = ys * sin(ecl/rad)

C  RA and Dec in degrees:
      RA = rad*atan2(ye,xe)
      Dec = rad*atan2(ze,sqrt(xe*xe + ye*ye))

      GMST0 = (Ls + 180.0)/15.0
      LST = mod(GMST0+UT+lon/15.0+48.0,24.0)    !LST in hours
      HA = 15.0*LST - RA                        !HA in degrees
      xx = cos(HA/rad)*cos(Dec/rad)
      yy = sin(HA/rad)*cos(Dec/rad)
      zz =             sin(Dec/rad)
      xhor = xx*sin(lat/rad) - zz*cos(lat/rad)
      yhor = yy
      zhor = xx*cos(lat/rad) + zz*sin(lat/rad)
      Az = mod(rad*atan2(yhor,xhor) + 180.0 + 360.0,360.0)
      El = rad*asin(zhor)
      day=d-1.5

      return
      end
